clear all;

# declare variables
filename = "llegadasregistro";
report = "~/Documents/ITBA/ss/proyecto_ss/plots/report_R";
clases = 7;
time_start = 8.0;
time_end = 13.0;

for i=1:clases
	cant(i) = 0;
	O(i) = 0;
	M(i) = 0;
endfor

#abro archivo para el report
report = fopen(report, "w");
fprintf(report, "Resultados para R\n\n");

# abro, leo el archivo y lo cierro
file = fopen(filename);
fscanf(file, "%c", 49);	#consumo el header
i = 1;
while (!feof(file))
	time_in(i++) = fscanf(file, "%f", 1);
endwhile
fclose(file);

n = i-1;

# computo los tiempos entre arribos de clientes
time_max = 0;
time_min = 1000;
for i=1:n-1
	times(i) = time_in(i+1)-time_in(i);
	if (times(i) > time_max)
		time_max = times(i);
	endif
	if (times(i) < time_min)
		time_min = times(i);
	endif
endfor
n = n-1;

# computo las clases
dif = time_max - time_min;
fprintf(report, "N\t\t%d\n", n);
fprintf(report, "Clases\t\t%d\n", clases);
fprintf(report, "Maximo\t\t%g\nMinimo\t\t%g\n", time_max, time_min);
for i=1:n
	for j=1:clases
		if (times(i) > ((j-1)/clases)*dif && times(i) < (j/clases)*dif)
			O(j) = O(j) + 1;
		endif
		M(j) = ((j-1)+0.5)*((time_max-time_min)/clases);
	endfor
endfor

for j=1:clases
printf("%d\t%f\t%f\t%f\n",j,O(j),M(j),O(j)*M(j));
endfor

# computo el estadistico de maxima verosimilitud
sumMaxVer = 0;
for j=1:clases
sumMaxVer = sumMaxVer + (M(j) * O(j));
endfor
lambda = sumMaxVer / n;
fprintf(report, "lambda\t\t%f\n", lambda);


# calculo la cantidad esperada para la exponencial
for j=1:clases
a = (time_max / clases) * (j-1);
b = (time_max / clases) * j;
E(j) = ((1-exp(-b/lambda)) - (1-exp(-a/lambda))) * n;
endfor



# TEST CHI CUADRADO
# H0: PROVIENE DE UNA EXPONENCIAL
SumChi = 0;
fprintf(report, "\nTEST CHI CUADRADO\n");
fprintf(report, "H0: El tiempo entre arribos de R proviene de una distribucion exponencial\n");
fprintf(report, "Clases\tO(i)\tE(i)\t\tO(i)-E(i)\t(O(i)-E(i))^2\t\t(O(i)-E(i))^2 / E(i)\n");
for j=1:clases
	DifCuad(j) = ((O(j)-E(j))^2);
	Chi(j) = DifCuad(j) / E(j);
	SumChi = SumChi + Chi(j);
fprintf(report, "%d\t%g\t%g\t\t%g\t\t%g\t\t\t%g\n",j,O(j),E(j),O(j)-E(j),DifCuad(j),Chi(j));
endfor
fprintf(report, "\nEl estadistico Chi Cuadrado es %g\n",SumChi);
fprintf(report, "Buscar en tabla el Chi Critico para %d grados de libertad y el nivel de confianza deseado\n", clases-1);
						  
# TEST Q-Q
y = sort(times);
for i=1:n
	gamma(i) = (i-0.5)/n;
	gammaInv(i) = -(1/lambda) * log(1-gamma(i));
endfor
#plot(gammaInv,y, '.');
#xlabel("Quantil Teorico");
#ylabel("Quantil Observado");
#print -deps "~/Documents/ITBA/ss/proyecto_ss/plots/PlotQQ_R.eps";

# histograma
hist(times, M);
#axis ([time_min, time_max, 0, 55], 'equal');
axis ("labelx", "tic");
xlabel("Tiempo entre arribos");
ylabel("Frecuencia");
#print -deps "~/Documents/ITBA/ss/proyecto_ss/plots/Histograma_R.eps";

fclose(report);